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ESTIMATION OF PARAMETERS IN ORDINARY 
DIFFERENTIAL EQUATIONS USING INTEGRAL FUNCTIONS 


by 


(C) YASEF OZSARFATI 
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IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE 
OF MASTER OF SCIENCE 


DEPARTMENT OF CHEMICAL ENGINEERING 
EDMONTON, ALBERTA 
PALL 19/73 


ABSTRACT 


Given a model which consists of a system of differen- 
tial equations describing a certain physical phenomenon, a 
study has been made of the problem of estimating the unknown 
parameters in the model from experimental data. In such 
problems, the integration of the differential equations 
describing the process and sensitivity coefficients is the 
most time consuming part. An attempt has been made to reduce 
the time requirements by introducing into the original model 
approximating functions of some kind. The resulting expres- 
sions do not require one-step or multi-step integration 
methods which are time consuming. 

It has been assumed that the only errors involved are 
those in the observations, and no dynamic disturbances are 
present in the system. 

Three problems have been studied; the first one is a 
boundary value problem and the last two are initial value 
problems. As mentioned earlier, the original model in each 
problem consists of a system of differential equations. Once 
this system has been transformed to a system of integral 
equations, functions approximating the observed variables 
in terms of the independent variable have been introduced. 
In the first two problems, the systems of differential equa- 


tions have been integrated using the Green's function approach. 
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The resulting expressions have been integrated analytically. 
In the third problem, which represents a more general case, 
where analytical integration is not possible, use of simple 
integration scheme, such as the trapezoidal rule has been 
made. 

The main problem associated with the approach employed, 
as will be seen, lies in the difficulties encountered in 
establishing proper confidence intervals for the parameter 


estimates. 
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CHAPTER I 


INTRODUCTION 


Finding a mathematical model which will describe a 
certain process, be it separation, chemical reaction, or 
a transport process, has been of growing concern in the 
field of chemical engineering. Sometimes the mechanism of 
the process is partially understood and an empirical rela- 
tionship is used to describe it. Sometimes complete under- 
Standing of the mechanism is possible, with the resultant 
mathematical model available. In both cases, complete 
information about the model is possible only by estimating 
the unknown parameters in it, based on observing the system 
behavior. The observations are usually made at discrete 
points, although in the case where time is the independent 
variable, one has additional access to continuous observa- 
tions. The estimation schemes to be reviewed below are 
non-sequential in nature, i.e., they apply to the case where 
estimation is done subsequent to observations. Particular 


references to kinetic studies will also be made. 


A. Literature survey 


Until recently, kinetic analysis of systems linearized 


by appropriate transformations have been common. These 
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approaches suffer two major drawbacks. First, they cannot 
represent real systems encountered in many situations, 
which are non-linear in nature, and second, under trans- 
formations, there is no guarantee that the original error 
distribution which has been judged appropriate will be pre- 
served. The availability of modern computers, on the other 
hand, has made the analysis of non-linear systems feasible. 
Examples for using analytical integration methods prior 
to estimation are given by Levenspiel (20), Wright (32), 
and Kittrell, Mezaki, and Watson (15). For simple systems 
One cannot object to this approach. However, it should be 
kept in mind that, although some systems will lend them- 
selves to analytical integration, the resulting models to 
be used for estimation are rather complicated to work with. 
The three-component, monomolecular, reversible reaction 
system analyzed in Chapter IV can be cited as an example. 
The general concept behind the estimation procedures 
is discussed by Bard and Lapidus (1), Quon (25), and Nieman, 
Fisher, and Seborg (23). They can be classified into three 
main groups; least-squares, maximum likelihood, and Bayesian 
estimation. The first, merely requires the assumption of zero 
expected value of the observation errors. The second, assumes 
a distribution for the observation errors, usually the norma] 
distribution. The weights, which express one's confidence in 
the observations of different variables, are considered to be 
part of the system and estimated together with the para- 


meters. The third method takes into account the existing 
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information about the parameters prior to the experiment. 
It is particularly useful for sequential experimentation. 

So far, only the existing criteria for estimation 
have been mentioned. An equally important aspect of esti- 
mation is the computational scheme employed to satisfy the 
criteria mentioned in the preceding paragraph. To this end, 


consider the initial value problem described by: 


Vou =s Ly yx say (I-1) 


y(X9) zi Yo (I-2) 


where y = observed dependent variable, 


x = independent variable, 


a = p-vector of unknown parameters. 

The prime stands for the derivative of y with respect to x. 
One could, as mentioned earlier, integrate Equation 

(I-1) analytically, if possible, and choose the resulting 


equation as a model: 


The alternative would be to work directly with the model 
given by Equation (I-1). In any case, the estimation pro- 
cedures to be outlined below are the same, regardless of the 
model being used. The only difference is that numerical 


integration might be necessary if one is working with the 
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model given by Equation (I-11). 

The existing computational methods can be divided 
into two main classes; direct search methods and gradient 
methods. Usually, direct search methods are used to provide 
an approximation as a starting point for the gradient methods 
which are iterative in nature. 

All the gradient methods provide a correction 6 to 
the parameter estimates at the iteration stage under ques- 
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where 2 = a scaler determining step length, 
3S ineP) ar at 
id = loa 3 Jo 28 SA | 
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S = objective function to be minimized, a function 
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Ee ay ps) 33 


y = N vector of observed values of y, 
y = N vector of predicted values of y by the algorithm, 


N = number of observation points. 
The different forms of the gradient method emerge as 
a result of the choice of A and R. This point has been dis- 
cussed by Bard and Lapidus (1) in detail. One of the most 
commonly used approaches is to express R in terms of the 


sensitivity coefficients which are the partial derivatives 


7 ewe Ts" 

behivib ed nao ebo ito a Ianotss tugn 
jnstbstp bas ebodten: do er devtha 
he ane 
2bodtom tnetbsie eft wet snfog. Perea? a ‘ote 


.Stutsa at ovisoye: 


sbiverq ot beay ax ebortem dow 


oF & noifasites 6 sbiverg 2zbodsom jasthere 94 


-29uUp "sdb Sspai2 notisyssf sdd its 2otnmttes 


T \ 2 
(fh~]) vH-; -— =| @ P 


tons qasz2z ontafrayeteb 191622 6 


c r < 26 
> re ee nd = 
l 
norjonvt s-,.besimintm 3d oF notsonut svitvostde = 


7 


: 
«3 2 To 


i> 


a | 2 at 
y 


-¥ to esulsv bevisedo to Yotsev H yom 


ivi 
oa 


ere ee, ae ‘ | ae 
;untttopi[s sd3i yd y to eeulev bstatbete to r10f29y wy : 
' | ie » 


.2intog notteviezde to vedmun = HW 


26 spvsme bodiom tnethave SAT Yo earr0t frorsttib odT 
-2tb need 260 dnteq erat .A bas & to et One onz to ste tS 
Jzom sant to on0 .Ttsteb at (1) eubtass rer bind vd bea 
sit to emist anf A 2eesvqxe oF 2f staat ln) el oem 


Visevineb Tetin8q one S78 aera setter naen- tana 
Neve ONT 7 oy 


hy : 


ie MV 7 
es, Os 


of the dependent variables with respect to the parameters. 
In the case where the model given by Equation (I-1) is used, 
this requires numerical integration. An alternative method 
is to evaluate them by differences, using slightly perturbed 
parameter values. It should be pointed out that the com- 
putational requirements are the same in both cases. 

Quasilinearization has received considerable attention 
Ge oe ra. “h/t.  linis method converts the ron-linear initia 
Or boundary value problem into a series of linear initial 
value problems. Donnelly and Quon (9) have reduced the number 
of equations to be integrated by considering the structure 
of the equation. 

Hwang and Seinfeld (12) have shown that a gradient 
technique similar to the one summarized by Equation (I-4) can 
be developed using the notions of quasilinearization. In 
this case, R = (alg) 7h, and v = ge} in specific, where G 
is the matrix of sensitivity coefficients. 

One of the shortcomings of the quasilinearization 
approach is that the convergence cannot be guaranteed. To 
assure convergence, Donnelly and Quon (9) have perturbed the 
observed data to generate a new set of data which will fall 
in the region between the observed data and that predicted 
by the algorithm. Seinfeld and Gavalas (26) recommend the 
use of polynomial fitting to the first several data points, 
or the use of the steepest descent algorithm for the initial 
stages of calculation. 


Invariant imbedding (8, 17) also appears as a powerful 
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tool in estimation. Optimal control concepts are applied 
in the minimization of the objective function over the 
whole trajectory. This method is sequential in nature, 
i.e., after each data point is taken, a new set of states 
and parameter estimates is obtained. 

The use of expansion functions in the solution of 
differential equations is quite common (6, 16). They have 
been also used in the reverse problem of identification 
(27,528)% 

The approach in this work will utilize expansion 
functions in estimating constant parameters in systems which 
have already been identified and can be represented by a 
system of ordinary differential equations. Cresswell and 
Quon (6) have proposed a method for attacking a problem which 
has also been considered in Chapter III. Lee (18) has used 
expansion functions in determining parameters which are func- 
tions of the independent variable x. 

Without a passing reference to the methods stressing 
the design of sequential experiments the review would not 
be complete. A number of papers have been published in this 
field (3, 4, 14, 22). It is generally recommended that the 
experiment be run under conditions where the determinant of 
eure is maximized, where G is the matrix of sensitivity co- 
efficients. The parameter estimates obtained from the last 
experiment are used in evaluating the new determinant of 
Gig over a grid of experimental conditions to choose the 


most favorable ones. Box and Lucas (4) have shown that under 
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certain assumptions this choice is equivalent to minimizing 


the volume of the joint confidence region of the parameters. 


By (Objectives of the Study 


The main drawbacks of the gradient methods and the 
quasilinearization and invariant imbedding algorithms outlined 
in Section A is the time requirements for the numerical inte- 
gration of the differential equations describing the system 
and sensitivity coefficients. 

As pointed out in the previous section, the use of 
analytical integration methods is limited, even if analytical 
integration is possible. 

The main effort in this work is spent in combining the 
features of the gradient and analytical integration methods 
in a way which will reduce the computation time requirements. 
To this end, approximating functions, depending on the nature 
of the problem at hand are used. In some cases, it iS neces- 
sary to specify these functions completely (i.e., estimate 
the coefficients involved) prior to the estimation of the 
parameters of original interest. In other cases it is 
possible to estimate both parameter sets simultaneously. 

The problems studied in Chapters III and IV enjoy the pro- 
perty that the system of differential equations can be put 
into a self-adjoint form. In cases where this does not prove 
applicable (Chapter V) a more general approach is taken. 

In this work, it is assumed that the model describing 


the system is already available. It is further assumed 
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that no dynamic noise is present in the system. Thus, the 


errors are limited to those in the observations. 
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CHAPTER II 


THEORY 


Consider the initial value problem! 
y Ves utxsajyebae= tty (x)ev(x23)sx) Tan) 
Yix es 9 | (11-2) 
aly Nieaie (11-3) 


where the prime denotes the derivative of y with respect to 
x, the independent variable; a is a p-vector of unknown 
parameters. The general solution of Equations (II-1) to 


Chi =) 1S) 


HX 7 {8),1y .X sa) (II-4) 


where 


For simplicity of notation, only the case with one 
dependent variable is considered. The above treatment can 
be generalized without any difficulty to handle cases where 
there are more than one dependent variable. 
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x 
p(x:a) = exp[f u(&:a)dé] (11-5) 
x 
0 
In Equation (II-4) the term z(v,n,x:a) accounts for the 
Tiieial “conditrons.iso thatywequations (1-2). and (11-3) will 
be satisfied. The procedure for determining F(x,&:a) is 
Given. i U3 I). 


Now consider the boundary value problem 


iy ae) ates a) yin (ex) re Safi yalixadies VEX Sia) 5m) (11-6) 
= yy (Xp) 7 p1¥(Xq) ay Or. (II-7) 
Moy (Xe) + Poy(Xe) = S (11-8) 


The general solution of Equations (II-6) to (II-8) is: 


X 
im 
yxy =f G(xses alpine). viasay.e)pCerayde 
X 
0 
( ) 
+ Bia a) r G(x,X,:a) + i ni = S G(x,x,:a) (II-9) 


where p(x:a) is again given by Equation (II-5). The proce- 
dure for determining G(x,&:a) is given in (31). 
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The general solution of 
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y(Xq) = vy (II-11) 


is given by 


X 
y(x) = p(x:a)7! f p(e:a)#(y(&)v(E:a) e)de + v (bt-12) 
XQ 
Now assume that a set of observations {Ya Xa 
i=1,2,...,N} is available. On the basis of those observations, 


the estimation of a is desired. Suppose one replaces y(x) 


by using an approximating function of some kind: 
y(x) = $(x:b) (11-13) 


where $(x:b) is the approximating function, and b is the q- 
vector of parameters in it. Replacing y(&) on the right-hand 
side of either Equation (II-4), (II-9), or (II-12), depend- 


ing on the problem at hand, by $o(&:b) gives: 


y(x) = p(x:a,b) (11-14) 


Replacing y(x) in Equation (II-14) by the observed values, 


One obtains the following error model: 
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far sbg suede ly 2m vase NK (II-15) 


Now the problem has been reduced to the estimation of a and 


b such that the least-squares objective function 
e4 (11-16) 


is minimized. It should be noted that, in some cases, due 
to the structure of the equations, the interaction between 
the components of a and b cannot be avoided, so that their 
simultaneous estimation is not possible. One could, then, 


define a new error model of the form: 


eC) hare pa ee, eg. N (1lb=-17) 


where c is a vector with a dimension smaller than ptq, whose 
components have been obtained by lumping the components of 
a and b. Depending on the situation at hand, c could still 
retain some of the elements of a. After estimating c using 
least-squares, application of a second least-squares proce- 
dure would give the desired estimates of a and b. 

Another approach would be to estimate b originally, 


using the following error model: 
Peres Sl eMes2t Re GON (11-18) 


Using b, the estimate of b, the following error model would 
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POR et etl pe gy tg N (11-19) 


If possible, the use of a two-stage estimation procedure 
should be avoided, since the estimates of a will somewhat be 
dependent on the estimates of b. 

Note that the approach discussed so far has been based 
on the assumptions that the differential equation can be con- 
verted into a self-adjoint form, and the analytical integral 


of 


exists. It has also been assumed that an analytical expres- 
sion for F(x,&:a) or G(x,&:a) can be found. When these 
assumptions hold, the above approach is particularly appeal- 
ing, since it lends itself to analytical integration. When 
these assumptions do not hold, one would have to resort to 

a more general approach. To this end, consider the initial 


value problem: 
y'(x) = hly(x),x:a) (11-20) 
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able. Integrating Equation (II-20) gives: 


y(x) = 


Now assumed that y(x) 


whose parameters have 


Replacing y(&) in the 
o(&:b) gives: 


y(x) = 


Dietes | Co es 


Ua eld CI engi OT) ung oy co (11-22) 
XQ 


is replaced by an approximating function 


already been estimated: 
y(x) = o(x:b) (iio) 


integrand of Equation (II-22) by 


X 


[o> 


RES COGls (11-24) 
XQ 


Replacing y(x) in Equation (II-24) by the observed values 


gives the following error model: 


x X as 
SL x MeeOex Drake anlehie(e.b ive a) ae (11-25) 


The above error model 


1S: CONTINUOUS IM x] hence, 1 can oe 


evaluated at as many points as desired. 


Note the basic difference between the conventional 


estimation methods and the above approach. In the former, 


the differential equations describing the system and sensiti- 


vity coefficients are 


integrated using one-step or multi-step 


formulas. The step size used in integration plays an import- 


: v 1s Ce 
*, Pan. 


tne at Nase at 
re if or se i, i 
o ny = af 3 = 
maaan 1s ap 
w a Game te i ee vie al ou : 
noftonut patismixotwges 16 +e batt at ae se 
-betsati29 need xbssy(s oven 2 


(S$S-1T) 


: or. m, ane a os a@e a et? a nll 
es- : = (x Ey. « 
(ES-I1) (g x)¢ ( N uae, annie | 
.™ 

on & > 


ud (SS-II) notssupa to Mareetae ads at (3)y ent: 


. ry 7 am ¥ 
7 7 


fe x 7 
(BS-~11) 3b{6:3.(d:3)6)d | tu = (x)y 
o* 
4 ay 3 wee 
29urey bavisedo ent yd (AS-I1) moftsupa ar (x)¥ entree 


‘a 


:I[shom +o1v9 patwol fot 7 Be seule 


Pe 
, 


‘ ; | 40a 
(38-11) s0U8 3 AR SPERM Dray 3 LBRO es hi 
| 0 
at 
‘i Me 


ad neo 3f ,s3ash 2x nt awountiaos 27 lobom 168t78. ovods oT | 
.bettesb 26 2tntog ynsm 26 ¥ bodeuts 
fsnorinsvnos sid asswied sans tett tb dtasd std stov : 
«YomIot st nl TOS CaNee euanis ott bons 2bortsm notis 
-t#tenee bas meteye sat eatdinsesb 2notisups Istinevettth ¢ 
qste~ft3lum yo qete-ano patev betsipeint S15 2tnstotttsea y 
-~Ivoqmi as 2eyseiq notisipstat ah bens sste is 
~ aa i A 


a” es 7 
hee ee YY ae 


15 


ant role in the stability of the solution of these differen- 
tial equations. Thus integration should be carried out using 
fairly small increments of x. With the new approach, on the 
other hand, the integrand can be evaluated at a certain 
number of x points throughout the interval, and a quadrature 
formula, or a simpler formula such as the trapezoidal rule, 
or Simpson's rule can be used for integration. The number 
of points at which the integrands are evaluated could thus 
be reduced considerably, since the integrals would not be 
as sensitive to the step size used. 

A brief outline of the non-linear least-squares proce- 


dure employed is given in the Appendix. 
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CHAPTER III 


A NONIDEAL FLOW CHEMICAL REACTION PROBLEM 


A. Statement of the Problem 
Consider the differential equation describing the com- 
position of a single reactant along the length of an isother- 


mal, nonideal flow chemical reactor: 


+26 ~ S ~ pa-g(c) = 0 , O<xs1 canes 
Pe} thx dx 
with 
4 | de oy 
c= TP + apn oa at x=0 CLT l= 2)) 
ce - Lf 
Ane Pid: @ Neier a (III-3) 
where c = dimensionless measure of the composition of 
the reactant, 
x . = dimensionless measure of the distance from 


the reactor inlet, 


Pe = Peclet number, vL/E, 

Da = Damkohler number, > k etge for homogeneous 
reactors, 

Q(c) = term accounting for the kinetics of the decom- 


position, assumed to be of the form cf, 


Vv = superficial velocity, 
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n = reaction order, 

L = tube length, 

E = axial dispersion coefficient, 
k = reaction rate constant, 


Co reactant composition at x=0"” 

The derivation of Equation (III-1) can be found in 
Levenspiel (21). The boundary conditions (III-2) and (III-3) 
have been discussed by Dankwerts (7), Wehner (29), and 
Pearson (24). 

A discontinuity at the left endpoint, Equation (III-2), 
is the mathematical approximation of a given physical situa- 
tVono@ Le 1s TOM tis reason that, Limits trom tne Loni mus Cc 
be considered near this end point. 

In Equations (III-1) and (III-2) Pe accounts for the 
nonideality in flow; Pe + 0 corresponds to complete mixing 
in the tube, Pe > » corresponds to plug-flow behavior. 

The aim of this chapter) 1s to estimate Da, Pe, and n, 
given a set of observations {C,.x,3 Tol, 2a. Gane where a 


j 


stands for the observed value of composition at x = Xa. 


Ba “Error Model 
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where & is a dummy variable introduced in order to carry out 
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Replacing Q(c) in Equation (III-4) by c" gives: 
een | Pe(x-&)_n 

edge I= Dabiec (ejde,+ fre Ce) del (I1I-5) 
0 X 


Let c"(x) be approximated by a polynomial of degree M: 
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c"(x) a_ x (T11-6) 
m 

Combining Equations (III-5) and (III-6), and replacing the 

left-hand side by Ca the observed value of composition at 


XoF Kae gives the following error model in discretized form: 
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The analytical integration of the terms in Equation 
(III-8) is quite straightforward. With the error model 
given by Equation (III-8), the objective function to be 


minimized is: 
N 2 
S, 5 ) Wie; (III-10) 


where Wes 1) ,c,.cagN are welgnts, itntroduced, to reflect 
the precision of each individual observation. 

Although the original purpose was the estimation of 
Da, Pe, and n, given a set of observations, the problem has 
now been transformed into one involving the estimation of 
Pe and bi m=0,1,...,M. In the related model described by 
Equation (III-8) Da and n do not appear explicitly; they have 
been lumped into the bi Ss This makes it necessary to use 
a two-phase estimation procedure consisting of: 

Phase 1: The simultaneous estimation of: 


(a) the effect of the nonlinearity in flow 
(accounted for by Pe), and 


(b) the effect of the decomposition rate 
Marccounted for “bY bm; m=0,l,..~2sM) on 
the composition profile. 
Phase 2: . The subsequent estimation of Da and n. 
| Phase 1 involves minimization of the objective function 
S;, given by Equation Gil-9).98 This, Wsesacconplatished byus.- 
ing the nonlinear least-squares technique. At the end of 


Phase 1, estimates of Pe and b> m=0.,1,....M are avallabie; 


these will be denoted by Pe and bi: MeO eles sy Me 
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From the inspection of Equations (III-6) and (III-9) 
10 follows. that: 


A 
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The required error model for Phase 2 is obtained by combining 


Paueations (1ll=+8) “and CLIT-11)): 
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Note that the model given by Equation (III-12) is 
linear with respect to Da and n, when transformed into log- 
arithmic coordinates. This procedure reduces the computation 
time required for Phase 2 considerably, compared to that for 
Phase 1. However, it is realized that the error distribution 


has also been transformed. 
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C. Generation of Data 


1. Generation of Noise-Free Data 

For n=1, the analytical solution of Equations (III-1) 
to (III-3) was used to generate the noise-free observations 
at eleven equally spaced points along the length of the reac- 
tor, corresponding to various Da and Pe values. The numeri- 
cal solution for n=2 was taken from Creswell and Quon (6). 


For n=l, the analytical solution is given by 
an = ky exp[Pe(1+8)x./2] + ky exp[-Pe(B-1)x./2]; 


{ealatbae eee i (III-14) 


where K, = ee ee fle ts} 


(B+1)" exp(B:Pe) - (8-1) 


z 2(8+1) £ 
ko ; Cee eS (111-16) 
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c. noise-free observation of c at x = Xa. 


The various data sets used are coded in Table III-1. 
The noise-free observations corresponding to those data sets 


are shown in Tables III-2 to III-/7. 
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TABLE III-1 


CODING OF THE VARIOUS DATA SETS USED 


Data Set Da Pe n( known) n( unknown) 
A 1.500 6.000 1.000 i 
B 3.000 6.000 1.000 s 
C 1.500 50.00 1.000 is 
D 3.000 50.00 1.000 aS 
E 2.000 6.000 Po: 1.000 


F 2.000 6.000 -- 2.000 
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NOISE-FREE OBSERVATIONS 
DATA SET A 
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TABLE III-3 


NOISE-FREE OBSERVATIONS 


DATA SET B 
X co” (x) 
0.0 0.73205 
0.7 0.58772 
0.2 0.47186 
0.3 0.37886 
0.4 0.30427 
0.5 0.24451 
0.6 0.19683 
0.7 0.15922 
0.8 0.13056 
0.9 0.241708 
1.10 0.10324 
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TABLE III-4 


NOISE-FREE OBSERVATIONS 


DATA SET C 
X c” (x) 
0.0 0.97167 
0.1 0.83989 
0.2 0.72597 
0.3 0.62751 
0.4 0.54240 
0.5 0.46884 
0.6 0.40525 
0.7 0.35029 
0.8 0.30278 
0.9 0.26175 
1.0 0.23272 
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TABEE INH-5 


NOISE-FREE OBSERVATIONS 


DATA SET D 
Xx c* (x) 
0.0 0.94627 
0.1 0.71240 
OF..2 0.53634 
O}..3 0.40378 
0.4 0.30399 
@)..5 0)., BASE 
0.6 0.17230 
O../ O}.. 229731 
0.8 0.09765 
0.9 0)..07353 
Liga 0.05832 
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TABLE III-6 


NOISE-FREE OBSERVATIONS 
DATA SET E 
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TABLE* I TI-7 


NOISE-FREE OBSERVATIONS 


DALLA, SET F 

x c* (x) 
OO 0.83129 
O24 0.74012 
OFZ 0.66559 
U3 0.60372 
G4 Ores obey 
OS 0.50764 
0.6 0.47017 
Oa 0.43862 

ac 0.41305 

a) 0.39476 
je 0.3872/ 
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2. Generation of Noisy Data 


* 
Each solution, Cis corresponding to a certain Da, Pe, 
and n was contaminated with the addition of random numbers 


to generate twenty sets of noisy observations, using: 


Ci TOG eet lace cost ls dalecm acd  (I1To 18) 
where a, = true solution at x = Xa 

Ci; = noisy observation at x = x,, for the jth run, 

ae = random number with Gaussian distribution, having 


zero mean and a standard deviation of 0.03, 
i. Tle) = Fy(el0, 0.0009)”. 


Subroutines GAUSS and RANDU, supplied by the IBM 1130 


scientific subroutine Package (13) were used for this purpose. 


D. Treatment of Data 


1. Weights, and.Initial, Guesses.of Parameters 


All the computations were done using equal weights. 
This choice was based on the knowledge of the distribution 
of the observation errors, since they were all generated from 
the same parent distribution. The weighting matrix used was: 
1 
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=(NXN) ~ 52 =(NXN) (III-19) 


the identity matrix, 
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For the minimization of S| Gewasgset equal to 0.03. For 


the minimization of So» o2 was approximated by 


(III-20) 


where p is the number of parameters estimated in Phase 2. 

For the noise-free observations, second, third, and 
fourth degree approximating polynomials were tried. For the 
noisy observations, attempts to use polynomials with degrees 
higher than two resulted in parameter estimates with high 
variances. The initial guesses of Pe and bi: m=0,1,...,M 
for the noise-free and noisy observations are shown in Tables 


Ili-8 and III-9. 


2. Treatment of the Results Obtained from Noisy 
Measurements, 


It was noted before that for every set of Da, Pe, and 
n, twenty sets of noisy observations were generated; the 
results of the individual runs were averaged and statistical 
conclusions were drawn regarding the means and variances of 
the parameter estimates. 

The nomenclature given below will be used throughout 


the remainder of this work: 
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p= 2! (111-23) 
J 
Pus! = least-squares estimate of parameter p, 
J obtained from the jth run, 
53 te least-squares estimate of the variance 


of p, obtained from the jth run. 


In words, is the variance of parameter p, obtained 


from the individual parameter estimates, whereas S° is the 


Teh IRs) 


variance of p, obtained by averaging the estimates of para- 
meter variances predicted by individual least-squares runs. 


p is the mean of the individual estimates of p. 


E~ Results and Discussion 


1. Noise-Free Observations 

In Tables III-10 and III-11 results for Da, Pe, and n 
are shown for noise-free observations. 

A study of these tables shows that the estimates of Da 
are closer to their true values for the case where n is con- 
sidered known than for the case where it is treated as an 
unknown. (The estimate of Pe would not be affected by this 
choice, since its estimation is carried out in Phase 1.) 
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TABLE III-8 


INITIAL GUESSES OF THE PARAMETERS 
NOISE-FREE OBSERVATIONS 


Degree of 
Approximating p fo) pio) p f°) p{0) p(?) 
Polynomial 
a 5G -5.0 5.0 -- -- 
3 5.0 -5.0 10.9 -5.0 -- 
4 5.0 -5.0 >.0 -5.0 520 
TABLE III-9 
INITIAL GUESSES OF THE PARAMETERS 
NOISY OBSERVATIONS 
Data Set p{0) pio) pf) pe (9) 
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polynomial degree. This, however, does not mean that the 
degree of the approximating polynomial can be increased 
indefinitely; the minimum variance criterion to be discussed 
later in this chapter would be the dictating factor in this 


choice. 


¢. Noisy Observations 


The results of individual runs from one to twenty for 
different parameter sets are shown in Tables III-13 to III-18. 
All the results were obtained using second degree approximat- 
ing polynomials. The summarized statistical results are 
given in Table III-12, in accordance with the nomenclature 
given in Section D.2eof thisecnapter. 

Note that in Table III-12, some of the results for Pe 
have not been shown for Data Sets C and D. Studying Tables 
IITI-15 and III-16 shows that the variance estimates for Pe 
fluctuate widely. In Table III-15 the estimates of Pe also 
fluctuate widely. Thus, any statistical conclusions drawn 
from these results would not have much significance. 

Studying Table III-12, the following conclusions are 
drawn for the range of investigation: 

(a) The estimates of Pe are highly sensitive to the 
observation errors. 

(b) The normalized standard deviation of the estimates 
Gbnld. Ree ho Ba, is relatively insensitive to the value 
of Da, but inversely proportional to the: value) of Pe. 


(c) The normalized standard deviation of the estimates 
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TABLE III-12 


STATISTICAL RESULTS FROM NOISY OBSERVATIONS 


Data Set A B C D E F 
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TABLE III-13 


INDIVIDUAL ESTIMATES 


DATA SET A 
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TABLE III-14 


INDIVIDUAL ESTIMATES 


DATA SET B 
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INDIVIDUAL ESTIMATES 
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TABLE III-17 


INDIVIDUAL ESTIMATES 
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TABLE III-18 


INDIVIDUAL ESTIMATES 
DAV SET F 
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2 
Pe? 


Pe, and inversely proportional to the value of Da. 


of Pe, (S V/2 


/Pe, is directly proportional to the value of 

(d) The normalized standard deviation of all three 
parameter estimates are directly proportional to the value 
Of Dx 

(e) (S624 °aDa is higher for the case where n is 
unknown than for the case where n is known. 

The last observation suggests that the estimate of the 
effect of the combined reaction parameters (i.e., Da-Q(c)) 
would be more reliable than the individual estimate of Da 
and n themselves. This limitation is not surprising when 
one keeps in mind the set up of the estimation procedure in 
which a polynomial is introduced to approximate the combined 
effect of Da and n. 

Note also that the estimate of Pe does not depend on 
one's knowledge about n; regardless of whether n is known 
or unknown the same estimate of Pe will be obtained in both 
cases. This fact is also due to the way the estimation proce- 
dure has been set up, which consists of two phases. It also 
makes sense physically; it cannot be expected that the know- 
ledge about a kinetic term (i.e., n) should affect a disper- 
sion characteristic term (i.e., Pe), whereas it is natural 


for 10. to atftect another kinétic term, toe... Ba. 


F. Comments on the Error Model 
(a) Note that the right-hand side of Equation (III-12) 


is a continuous function of x. This suggests that the evalua- 
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tion of the model error is not restricted to the use of the 
original x values at which data was taken. Thus, the number 
of points used in Phase 2 could be increased indefinitely. 
In the limit, the objective function S5 would be expressed 
as an integral of the model error over the interval (0,1). 


Equations (III-12) and (III-13) in continuous form become: 


] 
alee a eet e224) 
Ke , 


x) dx (III-25) 


The use of the above model is restricted to the case 
where n is a known integer; in this case, a quadratic expres- 
sion in terms of Da (as the only unknown) would be obtained. 
The choice of the reasonable root would be subject. to inspec- 
tion. Note that the algebra involved would get more and more 
complicated with increasing n. 

Equations (III-24) and (III-25) could be used to 


express Da as an analytical function of Pe and bee m=O) to. k.c sll, 
Da = F(Pe,b) (III-26) 


The estimates of the expected values, variances, and covari- 
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ances of Pe and bi: m=0,1,...,;M would readily be available 
from Phase 1. The following relationships would give an 


estimate of the expected value and variance of Da: 


we M Z 2 
] Or heed Tonk 622 
E(Da) = F(Pe,b) + > i —r, § + = S (III-27) 
é m=0 ab“ by é yPe- Pe 
V(Da) = ; (ave Coe Gaon ence 
m=0 2m Dh Be Fs 
aF_\/3F 
+ 2(55-) (Sp—) cov(by »b,) 
0 ] 
OF oF 
+ ATURE Re Ee (III-28) 
where E(Da) = expected value of Da, 
V(Da) = variance of Da, 
se = variance of Pe 
Pe : 
os = variance of b_, 
= m 


cov(b;,Pe) = covariance of b; and Pe, 


j K? = covariance of b, and by): j#k. 


All the terms on the right-hand sides of Equations (III-27) 
and (III-28) should be evaluated at the expected values of 
Perand sbi. .thatoisy Pe and b. The second and third terms on 


the right-hand side of Equation (III-27) could preferentially 
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be omitted without any serious loss of accuracy. Thus, only 
the calculation of the first derivatives would be necessary. 
(b) One could also consider solving the following 


equation for Da: 


] 
{ w(x)e(x)dx = 0 (III-29) 
0 


Although the computations involved in solving Equation 
(III-29) for Da would be much simpler compared to the com- 
putations involved in solving Equation (III-25), it should 
be kept in mind that the former sets the integral of the 
weighted model errors equal to zero, whereas the latter 
minimizes the integral of the weighted squares of the model 


errors, undoubtedly, a sounder criterion. 


G. Comparison with a New Model 


In order to test the effectiveness of the method, the 
results were compared with those obtained by estimating the 
coefficients of the approximating polynomial prior to the 
estimation of the parameters of interest, using the same old 


data. 


1. Error Model 


For convenience, Equation (III-5) is reproduced below: 


1 
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Suppose this time c(x), instead of c"(x), is approximated 


by a polynomial of degree K: 
. 1K es 
cexpaey }t ax ( Isl ies3)])) 
Combining Equations (III-30) and (III-31), and replacing the 


left-hand side by oF the observed value of composition at 


Ke aK as gives the following error model in discretized 


form: 
ei 
= K 4 m,n 
Sipe siete: LeyeDatie \ a) are rae 

0 m=0 

] POUX. eco IRT on 

to4448 land Guia dads 
Xe m=0 
j 
teby2matign (III-32) 


Suppose that the coefficients of the approximating polynomial 


have been estimated beforehand; call these estimates a3 


m=0,1,...,K. Note that K now represents the optimum poly- 
nomial degree, the estimation of which will be described in 


the next section. Replacing the unknown coefficients in 


ej; has been used in denoting both model errors defined 
by Equations (III-8) and (III-32), although they are different; 
this has been done for notational simplicity. Throughout the 
remainder of this work, the e's and S's should be interpreted 
only in the context of the model by which they are defined. 
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Equation (III-32) by their estimates, the error model 


becomes: 


# j K 
e: = Cc. - net Dalr{ the } fa Geet yn" dé 
0 m=0 
| Pelx.-&) KA 
Hoge dec GN ja, eyeeens 
Xie m=o ™ 
i . 
i Rd Pe cee tes | (III-33) 


The above model can be used to estimate Da and Pe 
simultaneously provided n is a known integer, since, other- 
wise the analytical integration of the terms in Equation 
(III-33) is not possible. Those terms will progressively 
get more complicated with increasing n. Note also that C, 
could have been replaced by its approximating polynomial. 

The above model was tested for n= 1. In this case, 


Equation (III-33) simplifies to: 
Xo. 1 
armedcfo-alpelpaomyalaipfase™de + f Ee eden 
i m 
0 0 Xs 


fal Zine ON (III-34) 
Again, unit weights were used in the minimization of: 
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2. Treatment of Data 

Forsythe's method (10) was used in generating poly- 
nomials in x to approximate the composition measurements. | 
The best-fit polynomial for the jth data set was chosen after 
observing the variance sae for different k values, ranging 


from two to five: 


N ‘0 . 9 
2 ee oi 
Si. = (Til 36) 
oJ N - (k+1) 
where Ci = observation at x = Xa TOY one. one run. 


Note that the method itself is based on the approxi- 
mation of a set of measurements by: 


k 
) Gra x) (7) 


where $,(x); m=0,1,...,k are functions of x. Using the 
appropriate relationships, the above expression can be trans- 
formed to a polynomial in x: 
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ed vale eapey| 


The motivation behind using Forsythe's method stems from the 
fact that ill-conditioning encountered with high polynomial 
degrees can be avoided due to the orthogonality of $m(x); 
m=0,1,...,k. Once they and their coefficients have been 
Calculated, the transtormaciron to (11) 7s trivial. 
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Cf he = corresponding value obtained by evaluating 
the jth approximating polynomial of degree 
K Sat ex f= x 2 
i 

A flexible variance criterion was defined; for the jth 
run: the polynomial with minimum variance , : was chosen to 
represent the data. In cases where ‘ ; kept decreasing 
slowly with increasing k, that value of k for which the next 
higher order polynomial caused less than about 10% reduction 
‘ Z ; 
in Sk oj was picked. 


The weighting matrix used is again given by Equation 


(III-19) with o = 0.03. 


3. Results and Discussion 


(a) Noise-Free Observations 
The results for noise-free observations are given in 
Table III-19. Again, note the improvement in the accuracy 


of the estimates with increasing polynomial degree. 


(b) Noisy Observations 


With data sets C and D, the original model seems to be 
superior as far assahe estimates of Re ‘ame jconcerned. The 
dependence of the estimates of Pe on the coefficients of the 
approximating polynomial cannot be overlooked, as was shown 
by a few runs with the same data set, but polynomials of 
different degrees. The results for Da, on the other hand, 
were quite good; even for ill-conditioned cases where Pe 
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was fairly close to the true value. This suggests that in 
cases where an accurate estimate of Pe is not available, 
the assumption of plug-flow behavior in the system will 
still give a fairly reliable estimate of Da. 

The summarized statistical results for data sets A and 
B are shown in Table III-20. The individual results are 


shown in Tables III-21 and III-22. 


H. Time Requirements 


All the calculations for this and the next two pro- 
blems were done using the IBM 1800 computer, with double 
precision. The time needed for computation was calculated 
by taking the difference between the time when computation 
ended and execution started. This included the time for 
building the coreload, reading the data, and printing. 
These, however, did not have an appreciable effect on the 
computing time, compared to the length of time spent on 
estimation. 

(a) Approximating polynomial with unknown coef fi- 
cients: average time required > 4 seconds/iteration 

average number of iterations: 6-7 (per run) 

(b) Approximating Meo mania with known coefficients: 
As mentioned earlier, coefficients of the approximating 
polynomials were estimated prior to the estimation of the 
parameters of interest. For this purpose, a program was 
used which calculated the coefficients and their variances, 
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TABLE III-20 


STATISTICAL RESULTS FOR NOISY OBSERVATIONS 
COMPARISON MODEL 


Data Set A B 
Da 1.509 3.004 
Pe Boas GL 578 
~2 
cae 0.0036 ~ 0.0154 
5° 2.082 1.813 
Be ; 
a 
ca 0.0010 0.0021 
52 W123 6.651 
o 
ae Wi 
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ua. 0.040 0.041 
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TABLE IITI-2] 


INDIVIDUAL ESTIMATES USING THE COMPARISON MODEL 


DATA SET A 
Run No. Pe Da Se so x10" 
] 8.437 B. S82 1.9996 1.026 
2 10.242 1.594 ROW fee 1.041 
3 5. 286 R542 0.6329 b.038 
4 4.629 §. 500 0.3379 1.005 
5 7.807 G.605 1.4700 1.054 
6 5.980 1.483 0.6925 0.977 
7 8. 32s 1.493 1.9734 0.994 
8 6.640 b..423 0.9750 0.939 
9 292 1.455 0.4882 Or. 915 7, 
10 5, R29 1.431 0.5856 0.966 
1] B Bike 1.574 1.4949 1.047 
12 6.790 bk. 482 O82 Q:., 91531 
13 4.982 v.43] 0.4290 0.958 
14 6.856 886 0.9845 1.039 
15 G., 2} 1.471 0.7518 0.978 
16 8.668 1.600 2 ROMS 1.070 
17 5. Shs7 1.500 0.4933 0.995 
18 6... 2901 13,500 0.7949 1.026 
19 6.029 2.457 0.7425 0% 983 
20 6.674 1.543 0.9630 E.038 
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TABLE III-22 


INDIVIDUAL ESTIMATES USING THE COMPARISON MODEL 


DATA SET B 
Run No. Pe Da San Sosa 
1 288 8.006 l S52 2.026 
2 ari9V'1 Saleli2 Terbr2e7 5 201013 
5 6.069 SeH 19 0.4909 21.268 
4 4 057/ 3.092 0.2679 25.5215 6 
° 1 shy) 311 9.0 0.8915 Zr 1516 
6 5.848 2.943 0.4256 2.016 
7 7.780 2.916 0.9480 1.948 
8 br: 953 2.787 0.4569 1.896 
9 6t.078 2.908 0s. 2986 2.04] 
10 5.439 2.860 0.3782 2.062 
11 7.869 Sy) iS 0.9396 2095 
1 6.218 2% 886 0.4974 1.884 
ie 4.834 2.884 0 121697 2.082 
14 6.886 32170 0.6328 2.160 
15 66:0 25 2.007 0.4696 2d Od 
16 nt ie StH 7? le2 77a 26.153 
17 54057 88 020) 0.3455 2.164 
18 6.239 86 02u 0h 5865 2nd 98 
19 ba 937 2.904 0.4712 2.071 
20 6.800 Bn. 09” 0.6661 2.215 
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Simultaneously. Time spent was about 8 seconds per run. 
For the estimation of Da and Pe: 
average time required : 4 seconds/iteration 


average number of iterations: 8 (per run). 


I. Qthers' Approaches 


Lee (19) has attacked the same problem, with n=2 (known), 
using a different approach; in this case Equation (III-1) 
isevbducea to two first order ordinary differential equations. 
Da and Pe are also treated as variables in x whose deriva- 
tives are equal to zero. The quasilinearization technique 
is applied to the resulting system of four ordinary differ- 
ential equations; this technique requires the simultaneous 
integration of four homogeneous differential equations with 
four different initial condition sets, to generate four homo- 
geneous solutions. Also a particular solution is generated 
by the integration of the system of four non-homogeneous 
differential equations. The initial conditions for the homo- 
geneous and particular solutions are chosen such that two of 
the four unknown constants appearing in the expression for 
the general solution vector coincide with Da and Pe. The 
other two can be expressed in terms of Da and Pe. Least- 
Squares procedure is carried out to find the values of Da 
and Pe which minimize a particular objective function which 
is non-linear in Pe. To this end, a random search technique 
has been used. Note that in Lee's approach non-linearity 


with respect to c(x) is dealt with by the use of quasilinea- 
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rization. Non-linearity with respect to the parameters 
(i.e., Pe) is taken care of by the use of a random search 
technique. In the present work, however, the only non- 
linearity to be considered is the non-linearity with respect 
to the parameters. Note, also, that the cumbersome inte- 
grations have been avoided by the use of a polynomial. 
Cresswell and Quon's approach (6) has been to fit an 
orthogonal Forsythe polynomial to the composition measure- 
ments, and obtain analytical expressions for Da and Pe in 
terms of the coefficients of the polynomial using moment 
relations. It should be noted that the availability of the 
analytical expressions is restricted to the case where n is 
a known integer. In this case, numerical experimentation 
is possible by using a random number generating routine into 
which the estimates of the means and variances of the co- 
efficients of the Forsythe polynomials are fed; using the 
analytical expressions for Da and n, estimates of their 
means and variances are obtained. For n being a non-integer, 
analytical expressions for Da and n are not available. 
Thus, the procedure described above for statistical analysis 


cannot be used in this case. 
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CHAPTER IV 
A THREE-COMPONENT MONOMOLECULAR REVERSIBLE 


REACTION SYSTEM 


A. Statement of the Problem 


Consider the monomolecular reversible reaction system 


Ay xy 
Koy kK 3) 
Kin = 73 
kK 30 
B (Hea a ESS 


oe ee eee 


kK53 


represented by: 


Xy> Xo and X4 are dimensionless measures of the com- 


positions of reactants A, B, and C, respectively. Ki 43 1,Jj= 


1,2,3, i#j, are the individual reaction rate constants. The 


differential equations describing the system are: 


xy = ~ (kay + kgy)X] + Ky oxo + k7 3X3 (IV-1) 
X5 = Koy Xz - (kyo + k3o)%o t+ kyg3X3 (IV-2) 
X3 = k31%X + KgoX% - (ky 3 + kog)X3 (IV-3) 
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subject to the initial conditions: 


x, (0) = X10 ‘ (IV-4) 
x, (0) = Xo 0 , (IV-5) 
x4 (0) = X39 (IV-6) 


In the notation above, prime denotes the time derivative of 
the variable in question. Since the compositions have been 


norimaltzed,-1.e.-3 
X, rex, tt xX, = Tee "eta 2) 10 (IV-7) 
then it follows that: 


xX, + X 5 tex 2 SAY "tb (IV-8) 
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Thus any two of the Equations (IV-1) to (IV-3) are sufficient 


to represent the system. Eliminating x, in Equation (IV-1) 


2 
and (IV-3) using Equation (IV-7) gives: 


Sc Sony? hi oue haa ge gue Kyo) aut Kas (IV-9) 
eee thee kes ee (ky a * Kog + Kay) X3 + ko (IV-10) 
Let 
k k k 
= Pa Ky; cs = Kas — oe K. ; (IV-11) 
12 23 31 
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and for simplicitys 
SK (IV-12) 


Ky» Ky» and K, are equilibrium constants, whose values are 
assumed to be known for the purposes of this problem. With 
the help of Equations (IV-11) and (IV-12), Equations (IV-9) 
and (IV-10) take the following form: 


Kans -[(1+K,)k, + k3]x, + (K,k. - k,) x, (IV-13) 


X9 (ka - Kok5)x, - [(K jk, + (1+K,)k,]xq + Kjk, (IV-14) 
Subject to the initial conditions: 
(IV-15) 
x .(0) =a ae (IV-16) 


The aim in this chapter is to estimate ky> kK» and k,; 
given the measurements of xy and X4 at different values of 


te starting trom C07, 


B. Possible Error Models 
The solution of Equations (IV-13) and (IV-14) in dis- 
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-a,t, t. a,(E-t,) 
ae eS ~ ! e [k, + b,x3(&)]dés (IV-17) 
-a3t, ti a,(&-t.) 
X30 § = X3 9@ + ; e [K,k, + b,x, (é)]dé; (IV-18) 
12152, N 
where Xj and X 9 ‘ are the dimensionless measures of the 


compositions of the reactants A and C, respectively, at 


t = t.. The remaining terms are defined below: 


a, (1+K) )k, + k, 


(I1V-19) 
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Note that although the original set of differential 
Equations (IV-13) and (IV-14) are linear with respect to the 
parameters, this is no longer the case for the solutions 
given by Equations (IV-17) and (IV-18). 

One way to attack the problem would be to approximate 
xy and X by polynomials in t, the coefficients of which are 
to be estimated simultaneously, with the original parameters 


of interest. To this end, let 
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Te he 
Ce) eRe = } pita (1¥=20) 
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Combining Equations (IV-17) and (IV-18) with Equations (IV-20) 
and (IV-21) gives: 
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Rearranging the integrands in Equations (IV-22) and (IV-23) 
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t. 
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The integration of the above expressions is quite straight- 


forward. The objective function to be minimized is: 


] 1 
we Pal ei oe ee 
] 3 
where a4 and a are estimates of the variances of the xy and 
X3 measurements. 

Attempts to minimize (IV-26) for the simultaneous 
estimation of Ky» Kos k3 and Pas Ais M=O legac cede EO ee oa 
using the models defined by Equations (IV-24) and (IV-25) 
would prove unsuccessful, due to the interaction between 


various parameters. In the following section, some approaches 


which will circumvent this situation are described. 
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1. Nonlinear Polynomial Model I 


Suppose the observations x ‘ and x 


are approximated by polynomials in t, so that: 
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(IV-27) 


(IV-28) 


Comparing Equations (IV-20), (IV-21) and (IV-27), (IV-28), 


note that in the former L and M represent arbitrary poly- 


nomial degrees, whereas in the latter, they represent the 


optimum degrees chosen according to the minimum variance 


criterion described in Section G.2 of Chapter III. 


Note also 


that Pi, and qn: m=0,1,....,L3; n=0,1,...,M are the least-squares 


estimates of pe and qn: 


Combining Equations (IV-17) and (IV-18) with Equations 


(IV-27) and (IV-28) gives: 
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will be called "Nonlinear Polynomial Model 


nonlinearity with respect to parameters ky 
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The model defined by Equations (IV-29) and (IV-30) 


cas 


Linear Polynomial Model 
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and k,. 


Introducing Equations (IV-27) and (IV-28) into Equa- 


tions (IV-13) and (IV-14) gives the following discretized 


model after integration: 
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This model will be called “Linear Polynomial Model" 
due to its linearity with respect to the parameters of 


interest. 


3. Nonlinear Exponential Model I 


So far only polynomials have been considered as candi- 
dates for approximating the measurements of x, and X3. 
Inspection of the composition profiles could suggest that 


the measurements be approximated by exponential functions: 


ae (IV-35) 
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(IV-38) 


where x and X34 are the equilibrium compositions of react- 
+) 


1,e 
ants A and C respectively. With this choice Equations (IV-35) 


and (IV-36) become: 
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reduces to: 


iu 
+ b1(%3 9 - X39) ! e e dé (IV-43) 


After performing the necessary integrations and rearranging, 


the following error model is obtained: 
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Model I". The least-squares objective function should be 


minimized to estimate Ky: Kos K3> Vis and v. simultaneously. 


3 
The original choice of xy and X3 to be approximated 

by exponential functions followed the inspection of the com- 

position profile generated with certain initial conditions 

and k values which will be given later. With other numerical 

values, the resulting composition profiles could dictate 

that other reactant sets be approximated by exponential func- 

tions. in any case, it should be Kept in mind that this 

approach is strictly dependent on the experimental observa- 
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4. Linear Exponential Model 
Combining Equations (IV-39) and (IV-49) with Equation 


(III-13) gives, after integration: 
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Equation (IV-46) reduces to: 
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The above model represented by Equations (IV-47) and 
(IV-48) will be called "Linear Exponential Model" because 
of its linearity with respect to the parameters of interest. 
Again, ky > Ko» k 35 Vi> and Vz are to be estimated simultane- 
ously. 


5. Nonlinear Polynomial Model II 


The models discussed so far, fail to take into account 
the effect of all three parameters simultaneously. The 
error model for xy is expressed as a function of kK, and K35 
whereas the error model for X3 is expressed as a function of 


k, and k.. In addition to that, the decoupling of vy and V3 


2 3 


eee ees oas eee a 
aie lust oe Oat te a 
pact igih is iy | a 
8s OTT 
ee eer at » 
(TR-VI) , Lalita tae +, ee “* 


i ae 
od ihe 


it: 
re 
st aan 


yey ) 


® as 


.3,V- ey 
a eee re 
iv (1% - or e4 * ot 
rane > - 
= ve" 9 ane ae 
(SPrNT) | Me naa Sel Sh a geen hg eee? geen 


w (Ra) I 


bos (V8-VI) anofseup3 ya bajnseayqet [sbom svods ” 
g92us29d "“{sboM [stinsnogxd yeantl" beff{so od [ftw (82- 
.tasisint to 21stsms 168q sat oF tooq2e07 Attw ys tosontl at 


-onsifumte bestsmtti2zs sd oF sts ev bie ef ¥ gt hich’ af 
4) 


Jnvos26 ofnt sist oF [tet .46?% G2 beeeureTb aide oaT 
eit .yfevosnstlumbe 21eteme1sq setdt tis vo doen 


IT [obo fsimonyfod rsenttnom a | 
, 
a! 


1g4 bas -A to notsonut 6 26 Pepys ee px wot Fh ot 
- . 
to notsonut 6 26 besezergxs at vor febom oe 


eV bas pv to pnt iquossb bisa mt Vicia 


aa a i‘ 


ia 


in the "Nonlinear Exponential Model I" limits the usefulness 
of this model, since one would not expect to obtain a unique 
set of parameters which will minimize the least-squares 
objective function in the region of interest. 

In order to express the error models for xy and X3 
as functions of Ky > Kos and K3> the following approach is 
taken. 

For convenience, Equations (IV-13) and (IV-14) are 


reproduced below: 
xy oT yaneas7 bs b 1X3 + ky (IV-49) 
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Differentiating Equation (IV-49) with respect to t yields: 
xy = a,x] + b,x (IV-51) 
Replacing x. in Equation (IV-51) by Equation (IV-50) yields: 
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x1' + a,x) 2 b, (bx, S)a,X3¥+ Knko) = F,(x:k) (IV-52) 


X, + = b,(-a,x, + b,x. + ky) = F(x:k) (IV-53) 
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The initial conditions of (IV-52) and (IV-53) are: 
x, (0) = xX 
X3(0) = X3 9 
(IV-54) 
AG LO) Saka Gg Meaty page Osa ay tak 
X3(0) = X3 9 = b3xX1 Q - 43X39 + Koko 


The solutions of Equations (IV-52) and (IV-53) with 


initial conditions given by Equation (IV-54) are: 
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Note that in the above model, both Equation (IV-57) and 


(IV-58) contain ky> Ko» and k This model will be called 
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“Nonlinear Polynomial Model IJ". 


6. Nonlinear Exponential Model II 
Replacing x, and x. in F > (x:k) by Equations (IV-39) and 
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Note that in the above model, both Equation (IV-62) and 


(IV-63) contain k k k Vie and V3: This model will be 
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called "Nonlinear Exponential Model II". 

The approaches taken in developing the various error 
models above can be generalized to an n-component monomole- 


cular reversible reaction system; in this case, the error 


model would consist of n-1 equations. 
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C. Generation of Data 


1. Generation of Noise-Free Data 
Data was generated by using the fourth-order Runge- 


Kutta integration scheme, with initial conditions: 


X10 = 1.0 
Xo 0 = 0.0 
X30 = 0.0 


The k values were taken from Wei and Prater (30) who have 
considered the reaction between l1-butene, cis-2-butene, and 


trans-2-butene in specific: 


Koy = 10.344 ; Kio = 4,623 
K35 = S30 HOn ko3 = 3,371 
ki 3 SY ip lic OO Ole ees Kk 34 Sy Seo 


With the terminology given by Equation (IV-12), the true 
values of parameters of interest are those ones given on the 
second column above. 

Integration was continued until equilibrium conditions 
were met, although only the first twenty observations (exclud- 


ing the initial conditions) were picked for estimation pur- 
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poses. Those are given in Table IV-1. 
With the initial conditions and parameter values given 
above, the generated composition profile for xy is mono- 


tonically decreasing, whereas that for x, iS monotonically 


3 
increasing. The composition profile for X3 passes through 
a maximum before reaching equilibrium. This is the reason 
for choosing xy and X3 to be approximated by exponential 


functions, rather than Xo. 


ew) Generation of Notsy Data 


Twenty sets of noisy observations were obtained by 
the addition of random numbers with Gaussian distribution 


to the true solutions for xy and X. using: 


§ mg.) Meee S al Gea 
Roeae X3 i i: Wis bee sa co (IV-65) 
h * * os . 
where Ry 4% 3 4 = true observations at t = tis 


= noisy observations at t = t., for the 
Stn. CUMN, 


2€5 = random numbers with Gaussian distribution, 
J both having zero means and a standard 
deviatiobeot 0.01% 


The initial and equilibrium conditions were not sub- 
jected to noise. Thus the terms X19? aS ek X39? and X36 


encountered in the earlier sections should be interpreted 


novtp 2sulsv ear: 
~onom a rig Abe a | morsel 
¢ff62taosonom # pt Ply ne : ; 
douordds 292289 ex 10% ott on 


noesst sat 2f etadt aur 


yd bantstdo sew inctiaennae ‘Veron to aoe vinowT a : 
nottudivietb nsteevsd dttw evediiun mobnsy to notdtbbs. a 
rontep x bae -x 10% anottufoe urd si "7 
wt id Rin 


ia 
a 


a 
: 
( 


(S8-VI) | : TG + t. + es tte C 


i) 


bm * pa ae eh i 
(2a-VI) OS iss sae f=e.f Ps Tie + re” th.&* (oe i, 


* *_ it a 7 
1;t = 3 46 anoftsvisedo suyt « hes tei a 


ont voY ..3 = ¢ 36 enotisvyeedo ee = x 
, biel AIL the e *oLhll 


Nofiudinseth nstezeusd difwia 
bysbnst2 6 bas 2ensem oves 
10,0 to, 


-duz ton $19W bine: ‘muivdtt ups bas ora 
o.e% DNB 4g) gk oP Pe ee 
betsiqisint sd ot tae ots 


Pet i 


ee : 
aa : a ae Sr oe 


TABLE IV-1] 


NOISE-FREE OBSERVATIONS 
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as noise-free initial and equilibrium conditions. 


D. Treatment of Data 

Forsythe’s method was used to fit polynomials in t to 
all the observation sets. The best-fit polynomial for each 
set was chosen according to the variance criterion described 
in Section G.2 of Chapter III. When evaluated at t=0, those 


polynomials did not yield the initial conditions (x, qo a0 ; 


a a OPO} ceXaecly... in order co ensure: tliat tne ane al 


conditions were satisfied to an accuracy of five digits, a 


weight of 10,000 was given to both X79 


and X30 in the 


Forsythe fitting procedure. Unit weights were assigned to 


the remaining twenty points. 


The initial parameter guesses are shown in Table IV-2. 


The weighting matrix used was: 
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I the identity matrix. (I1V-69) 


ey, “a 
ceooii toni 
mo 


dos yo? rcwnetes sit-sent 
badi1329b noriat ty sanétisy 
aie he 
e203 Ot 36 botsulsve nodW . 127 ebhed> 10 a 
.O.f = 0. rx) enotafbaes ratetat ‘odd brane don “bi | 
lettint sid todd srxwene of +ab%0 nt U1390K9 hee: ~ 


Pim 7 iyi 


6 ,2ttptb svtt to yostwoos me oF bSttetise ‘ovew anotdt 
ads nt o.e* bas 0.1% Atod of ‘Havin 26W 000, or ves ao) 


ot banotees siew etdptew sin .syubesorg pntasny on: ‘a | 


bash 


.S-VI ofdsT at nwonde s16 292z2oup Ted omsr1s5G [sidint ont 4 


.2tnfoq wand ontnt 


:26W bs2eu xiyisw outadp fon Oh i 


8 (nxn) 
[Pray * (usxns)# 
| xn! 9 
\ 4 I a 
(\d-V1) sa yh a hy 
[ 
i 
g? 
@h : 
(0a-VI) Men ee ode! 


> CM a Ny aviyesars -_ 


TABLE IV-2 


INT TTAL,GUESSES,OF THEsPARAMELERS 


Polynomial Models Exponential Models 
seas 5.0 
530 5). 0 
a7) Sd) 
--- 10.50 
ae ae %.0 


80 


08 \ a 


¢ ¢ me | . Fy | 


2AITIMARAY BHT wo eaeeaya 2 sa 


ae gen een he he a 


me et ee te Nm Se a 


eTabom is finanoqx3 sfabor fotwonytoR 


> ea ep a ge 


O.¢ O.¢ 


0,0f be total 


8] 


For the noisy measurements, both 1 and Oz were set equal to 
0.01. With these values, Equation (IV-66) reduces to: 


LAG’ ol (IV-70) 


WC ONX2N) = Le onxeNn) 
In other words, equal weights were used in the minimization 
of the objective function, since it was known that both 
observations had the same error distribution. 

For the noise-free measurements oy and O 3 were set 


4 


equal to 1.0 x 10 © somewhat arbitrarily. 


E. Results and Discussion 


1. Noise-Free Observations 
For convenience, the true values of the parameters are 


reproduced below: 


ky = gO oy k, So See Re ena C4 

With the polynomial models, fifth, sixth, and seventh 
degree approximating polynomials for x1 and X2 were tried. 
The results are presented in Table IV-4. A study of these 
shows the effect of increasing the polynomial degree on the 
accuracies of the estimates. One should be cautious, though, 
in using polynomials with high degrees in the case of noisy 
observations. 


Attempts to estimate the parameters using the Nonlinear 


; ae 
ot “Taupe rey ot9W e? ® ves 
| 203 z99ubs7 oa 


| oe Aion 
in) Pr 1 
iecrite ge 
notisstmintm edd at Bosal | 
dtod sed nwond enw abn | 
notiudratath wore pany aay > 
t92°979W' od ba Ps atnomoiwensm sales 


UE Psd Fdvs aan ccull 


ye 


Nofeat Id bogs 


anotisvy92ed0 satan “im 
| 


S16 218Jems16q Sit To zoulsv sunt ont .s9netnewnoD 


:woted 


ir 


yy, 


Sv.€ = od. IKE. + od, BSBIE & a oe i, 


itnsvez bane ,Atxte ,Adtt? ,e2efobom Istmonylog sad astw in 


ae 


ay i 


| | 
(HQuOKS ,2votsusS od biuode sH0 i2etemtses oft wan ri 
nt oe 


Y2fon to sz69 oft nt agornpeb Nota atiw Nene tek 
ie oie, ; 


Park 


.bst1t 9v9w eX bn6 xX TO? ef simonyfog pntsemrxorggs 9 
s2odd to ybute A .h<VI stdsT nt betnsesgtq 946 ahd 


ait no seipsb [stmonyloqg edt patesevant to toatte 


vsentinow odt pntauv erestemeveq edt stemtses of P ‘. isdtA . 
pe nN | 


ij hi . 
- ny \ 


re . a ae ia) f 
ons - vie i ‘ ‘li pan 


82 


Exponential Model I failed. This could be atrributed to the 


fact that Equation (IV-44) contains only k k and vz, and 


i? 
Equation (IV-45) contains only kos k,, and Vi: 


ai 


The Linear Exponential Model converged to a unique 
set of parameter estimates, but these were somewhat far from 
the true values. 

Among the exponential models the most successful one 
was the Nonlinear Exponential Model II. In this model both 
Equation (IV-62) and Equation (IV-63) are functions of all 
the five parameters. Thus, the parameters are completely 
coupled. 

One can conclude that the exponential models are not 
as successful as the polynomial models. This can be attri- 
buted to the lower flexibility of an exponential function of 
the form given by Equation (IV-39) or (IV-40) in approxi- 
mating a given set of data. Once the initial and equilibrium 
conditions have been fixed, there is just one exponential 
parameter left to be estimated. In the case where a poly- 
nomial is used as an approximating function, the number of 
parameters to be determined is equal to the degree of the 
polynomial if the initial condition is fixed. 

The results for noise-free measurements with the 


exponential models are presented in Table IV-3. 
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TABLE IV-3 


PARAMETER ESTIMATES FROM NOISE-FREE OBSERVATIONS 


Model 
Linear Exponential Nonlinear Exponential II 
4.264 4.582 
8.306 3.094 
3a 99 3.842 
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2. Noisy Observations 


Considering the relatively inferior results obtained 
with the exponential models, especially the linear exponen- 
tial model, using the noise-free observations, it was felt 
that further testing of these models with noisy observations 
was unnecessary. For the polynomial models, the results of 
the individual runs from one to twenty are shown in Tables 
IV-6 and IV-8. The summarized statistical results are shown 
in Table IV-5, using the nomenclature given in Section D.2 
Sreciapcey a) I, 

Studying Table IV-5 the following conclusions are 
drawn: 

(a) The variance estimates of the parameters are 
nearly the same with the three models. Based on this obser- 
vation, the null hypothesis that the variance estimates of 
each parameter, obtained with different models come from popu- 
lations with the same variance was tested. The results, 
obtained using an F-test for the homogeneity of variance, 
showed that this hypothesis can be accepted quite reasonably. 
From the nosie-free measurements, it is also known that the 
mean values of each parameter, obtained with different models, 
come from populations with the same mean. The above two 
points suggest that it is quite reasonable to assume that 
the parameter estimates obtained with the three different 
models come from the same population. Based on this 
assumption, the use of the linear polynomial model is 


strongly recommended, due to the simplicity of computation 
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STATISTICAL RESULTS FROM NOISY OBSERVATIONS 
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TABLE IV-6 


INDIVIDUAL ESTIMATES 
NONLINEAR POLYNOMIAL MODEL I 


ke k, k oh x10 
801 3.801 3.419 0.375 
326 3.180 3.969 0.300 
295 2.897 4.124 0.289 
548 3,397 3.825 0.349 
562 35262 3.969 0.315 
085 4.101 2.928 0.464 
533 3.077 3.022 0,363 
459 3,025 3.985 0,828 
624 3.88] 3.720 0,353 
012 3.583 3.540 0, 387 
624 2.818 4.247 0.296 
411 2.830 3.975 0.300 
64 2 3.348 3.042 0,330 
489 3.870 3.808 0.349 
BAN 4.744 2.619 0. 559 
183 2.993 4.206 0.276 
177 2.690 4.432 0.270 
508 BAO" 3.943 0,./31'8 
721 3.417 3.495 0,. Biby 
89] Bimoe 3..f550 0.388 
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TABLE IV-7 


INDIVIDUAL ESTIMATES 
LINEAR POLYNOMIAL MODEL 


Run No. k, k, k Si, X10 3 x10 sf x10 
2 3 

1 4.711 3...675 3.520 0.109 0.432 O.SY 
2 4.273 3%'028 4.087 0.088 04 388 04707 
3 4.236 2.001 4.157 0.087 0.368 Oe 207 
4 4.459 3.291 SSaedo ke: On. 107 0.446 US hee 
5 4.545 3.04] 4.15] 0.089 OvoZ7 0.254 
6 4.894 4.073 3.019 039129 0.501 Ay 5 
7 4.423 35427 3.615 0.101 0.430 0.298 
8 4.408 35097 3.998 0.099 0.396 0.289 
9 676 3.463 38.6258 Or1T4 0.444 0.324 
10 905 3.465 3.677 Ort 07 0.355 02296 
1] 4.310 22627 4.341 0.090 01332 Ox 256 
12 4.444 2.983 3.789 0.4.08 0.375 0.295 
is 4.448 32294 3.643 0.100 0.400 Os29" 
14 4.412 3.475 3.749 0.106 0.479 87313 
16 5.080 4.628 22802 0.146 0.618 0.434 
16 4.161 2:859 4.308 0.081 0.364 0.240 
ae. 4.068 21515 4.559 0.081 0:359 04299 
18 4.414 2.981 3.993 0: 095 0.350 0.266 
19 607 3.296 34 507 0.106 0.389 304 
20 4.664 32280 3.809 0.106 0: 390 0.304 
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TABLE IV-8— 


INDIVIDUAL ESTIMATES 
NONLINEAR POLYNOMIAL MODEL II 
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involved, since no iteration is required. 

(b) The sensitivity of the parameters to observation 
9 to k3 and kK). Dt 
should be noted that this observation would not generally 


errors goes in decreasing order from k 


be valid if initial conditions different than the ones in 
this analysis were used. As pointed out by Wei and Prater 
(30) and Donnelly and Quon (9), different initial conditions 
would lead to different parameter sensitivities, and certain 
initial conditions would not lead to any estimates at all. 
Therefore, the dependence of the parameter estimates to the 
experimental conditions should always be investigated in order 
to obtain a reasonable degree of reliability. 

(c) The actual variances of the parameter estimates, 


55 ; i=1,2,3, are larger than the averages of the individual 


cee estimates, 551% imld2t3 wiThiseiseattesulds Of the 
fact that the ee eee of the approximating polynomials 
are not considered as random numbers during the estimation 

of the parameters. Thus the confidence obtained in the para- 
meter estimates appears to be higher than it actually should 
be. A suggestion to obtain more realistic confidence inter- 
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F. .Time Requirements 


1. Nonlinear Polynomial Model I 


Average time required » 2 seconds / 1 teration.. 


Average number of iterations : Super Sane. 
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2. Linear Polynomial Model 


Computing time was very short in this case due to the 


linearity of the model with respect to the parameters. 


3. Nonlinear Polynomial Model II 


Average time required wo21 seconds/aterations 
Average number of iterations : 5 (per run). 
The time required for the estimation of the coefficients 


of the approximating polynomial was about 10 seconds per run. 


G. Comments on the Polynomial Models 


(a) Both the nonlinear and linear polynomial models 


59 and k, are discrete error 


models. If x, and X 3 are replaced by their approximating 


used in the estimation of Ky» k 


polynomials, the resulting models will be continuous in t. 
With this choice, the linear polynomial model takes the 


following form: 
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After rearranging, the above equations can be written as: 


L a re 
's N aN m 
e,(t) = X10 k,t + ay (1 + ra t)p,t 
M oq | 
te egret (1V-73) 
] "20 m+ 
Lp 
m m+ ] 
e3(t) = - X3 9 - Kokot - b, a rae 
M a. ss 
, Lo Ob Gp anet- “dat Sheree 
m= 


With this continuous error model, the least-squares objective 


function is expressed as: 
2 
€,(t)]dt (IV-75) 


A natural choice for te would be ty. i.e., the time at which 
the last measurements of X, and X3 were taken. 
Now one can set the derivative of S with respect to 


k k and kK. equal to zero and can solve the resulting 


Tate. 

expressions for these parameters only as functions of the 
coefficients of the approximating polynomials. An estimate 
of the expected values and variances of these coefficients 


would readily be available from the original fitting proce- 


dure. One can then use relationships similar to those given 
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bysEquationsy (112-22) ande(Libha28) in’ SectioneR of Chaptén 


III to obtain confidence intervals for the parameters, which 


would be more reliable than those given by the linear poly- 


nomial model itself. 


(b) The discussion would not be complete without a 


passing reference to the moment method. 


For convenience, 


Equations (IV-13) and (IV-14) are reproduced below: 


(IV-76) 


(IV-77) 


Replacing xy and X4 by their approximating polynomials, the 


following error model is obtained: 


Since the estimation of three parameters 


moment conditions require that: 
aus uly ee ar 
S [ 7 Ee, (t) ae e,(t)]t dt = 0; 
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With this method too, the sensitivity of 


is desired, the 


j=0,1,2 (IV-80) 
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observation noise can be estimated as explained in Part (a). 
Note that by expressing the error model in continuous 
form, one puts his whole confidence on the expected values 
of the polynomial coefficients, and does not worry about the 
experimental observations any longer. Thus, one might con- 
sider the discrete error models which are expressed in terms 
of experimental observations to be more realistic. The 
advantage of the continuous models is that they lend them- 
selves to subsequent statistical analysis concerning the 
estimates of the parameters. It should be noted that this 
is restricted to systems of linear differential equations, 
or nonlinear differential equations which can be expressed 


in self-adjoint form, as encountered in Chapter III. 


H.. Others’ “Approachés 


Wei and Prater (30) have studied the problem considered 
in this chapter and a family of related problems in consider- 
able depth, and have come up with a method which makes use 
of the "straight-line reaction path". They construct examples 
in which the set of parameters estimated from a single reac- 
tion curve will reproduce the original data adequately, but 
will give quite different curves than those expected when 
different initial conditions are used. Their suggestion for 
overcoming this inconsistency is either to run the original 
experiments with widely differing initial conditions so that 
most of the reaction simplex will be covered, or to locate 


the straight-line reaction paths, which, they consider easier. 
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It should be noted that this approach also necessitates the 
running of more than one experiment to locate the straight- 
line reaction path; once this has been done, no use is made 
Ofeene data: obtained prior: to that stage... In 'short, “their 
approach is one of a design of an experiment, rather than 
being an approach to the treatment of a specific set of data. 
In their work, attention has also been drawn to the 
hazards of using curve-fitting techniques for cases in which 
more than one exponential term is involved in the expression 
for the general solution. The model criticized is one 
obtained from the general solution to a set of linear differ- 
ential equations, which, in the case of a three component 


System reduces to: 
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Assuming that the initial conditions are known exactly, as 


has been done before, the following relationships hold: 


(I1V-83) 
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Note the close resemblance between this model and the 
Linear Exponential Model described by Equations (IV-47) and 
(IV-48). 

Wei and Prater have also given the analytical solution 
for the three-component reaction system described by Equa- 
Dronse (lv= 1) “to SUGV-6 )easwa tUnC tlONMO tek Som whu VLeWulOrE 
that fact, one might be tempted to use this analytical model 
for the purpose of estimating the parameters. Examination 
of the solution, though, leads one to conclude that it depends 
on the parameters in a rather complicated manner; the possibi- 
laty of using this .anatytical model is, thus, discarded. 

Donnelly and Quon (9) have attacked the same problem 
uSing quasilinearization. They have used different combina- 
tions of data obtained with seven different initial conditions. 
Naturally, the best results were obtained when all the seven 


data sets were used simultaneously. They report that using 
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only one data set generated with the initial conditions used 


in this chapter leads to acceptable results. 
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CHAPTER V 


A NON-ISOTHERMAL DECOMPOSITION PROBLEM 


A. Statement of the Problem 


The pyrolysis of ethane in the temperature range 1200 


to 1700°F is represented by: 
Gr hiper 2 Ce Heater (V-1) 


The reaction is carried out in a steel tube with no internal 
obstructions (e.g., catalyst), contained in a furnace. In 
the development of the relevant differential equations des- 
cribing the composition and temperature profiles along the 
length of the reactor, the following assumptions will be 
made: 

(a) Any pressure drop along the length of the tube 
can be ignored, 

(b) There is complete mixing in the radial direction 
and no mrxing in thevaxial direction (1.¢6., the flow is 
governed by plug-flow behavior), 

(c) No temperature gradient exists in the radial 
Girection, 

(d) Ideal gas law holds. 


With the above assumptions, the differential equations 
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obtained by setting up material and energy balance equations 


over a differential volume element of the reactor are (5): 


dz SO0C ATP kiln (v-2) 
dx ngR Tp +2 ij 
LER dz 
dT n AHR x 
R y (V-3) 
dx 1-z)c TieZCe oot 
Pa Re 10P.¢ 
where A' = cross sectional area of tube, sq. ft., 
Cy = specific heat of CoHe CoH, and H, correspond- 
ing to the subscripts A, B, and C, BIU/1bmole-°R, 
k = reaction rate constant. sea We 
aioe feed rate of ethane, lbmoles/hr, 
P = mean pressure inside the reactor, psia, 
q = heat input from furnace, BIU/hr-ft. of tube, 
R = jdeal gas constant, 10.73 psia-cu. ft./lbmole-°R, 
Tp = measure of absolute temperature, °R, 
Xx = measure of the distance from the reactor inlet, 
pte, 
Z = measure of the fraction of ethane converted to 
ethylene and hydrogen, 
AHp = heat of reaction, BTU/Ibmole. 


The purpose of the original problem given in (5) is to 


generate the conversion and temperature profiles along the 
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reactor, given the following information: 


meg 10 
ri 
ko aeoegod x 101°e SPs 1 (V-4) 
and 
Cy cal/gmole-°K 
: -3 -6.2 
CHE So.) Deets Ot x 10 Ty a ehO).: Liga 1:0 Ty 
CEn RM See Su ePORLoMRe OF oT. 2\'6. 88% 110 To (V-5) 
2 4 ° ° K ° K 
‘ -5 -6,2 
H, 7 POO Peat gs nomex 0 Ty =. KO) tee Ty 
i -3,+2 2 
AHp oT Oo Cat 8.50(T,-298) - 5.942 x 10 (T,-298 ) 
-6,+3 5 
+ 1.28 x 10 (T,-298 )J, BTU/1bmole (V-6) 
where Ty is the measure of the absolute temperature in °K. 
As seen from (V-4), k in general can be expressed as: 
we 
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Combining Eqwations (V-2), (V-7), (V-8), and (V-9) gives: 
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where D = ea ie Xx 10! (V-11) 
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In the derivation of Equation (V-10), use has also been made 


of the following relationship: 


ul = eS (V-12) 


The initial conditions for the differential equations des- 


cribed by (V-10) and (V-3) are: 


(V-13) 
and 


Tessar eT ;S. atm ox=0°: | (V-14) 


The estimation problem can now be posed as one of esti- 
mating a and ay given the thermodynamic information supplied 
by (V-5) and (V-6), together with the conversion and tempera- 
ture measurements along the length of the reactor. In other 
words, given the thermodynamic information, the estimation 


of the parameters related to kinetic effects is desired. 
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B. Error Model 
The general approach in the previous problems was to 
introduce some approximating function into the differential 
equations and try to estimate both the coefficients involved 
in the approximating functions and the original parameters 
Simultaneously. Specifically, when ordinary polynomials in 
xX are considered as approximating functions, the approach 
mentioned above will not be practical due to the resulting 
expressions for which the analytical integrals do not exist. 
This is the most general case encountered in the estimation 
of parameters from systems of ordinary differential equations. 
Below, a new approach to deal with this case will be described. 
Suppose the noisy conversion and temperature measure- 


ments have been smoothed by: 
m 
Eka = ag aes (V¥-15) 
and 
M m 
eG eh, Un iaae orn (V-16) 
respectively. L and M are the optimum polynomial degrees 
chosen in accordance with the minimum variance criterion 


described previously. 


Now, integration of Equations (V-10) and (V-13) gives: 
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xX. 
rn j 
pe Da, ; f(T,z:a,)dx; (V-17) 


ri x, a - AHpDa,f(T,z:a,) 
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where Tp has been replaced by T for notational simplicity. 
z. and le are the observed values of z and T at x = x,. 
Substitution of Equations (V-15) and (V-16) into Equa- 


Pons Uved/) “and (V=18) yields: 
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In interpreting Equation (V-20) it should be kept in mind 


that AHp, ¢, » C ; and ¢ are now functions of R(x), al- 
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notational simplicity. 


The least-squares objective function to be minimized is: 


e 


1 2 
a cael 
L- 


eT i (V-21) 


Lt 
of 
; 
2 é A : 
where 5 and OF are estimates of the variances of the conver- 
sion and temperature measurements. 


Note that the integrations required for evaluating the 


J€ 
model errors and the sensitivity coefficients CET eae. 
j ‘ 


J€ , 
i Tol 25. capone sg Cannot DeEReVa luatedeanalytically. 


aa; ex 
Instead, the integrands can be evaluated at a certain number 
of points, with the subsequent application of a quadrature 
formula or a simple formula such as Simpson's rule or the 
trapezoidal rule. The application of a quadrature formula 
usually requires the evaluation of the integrands at points 
other than those at which data was taken. In order to apply 
such a formula, z and T should be replaced by P(x) and R(x), 
so that the error model and the sensitivity coefficients are 
now continuous in x. A much simpler procedure would be to 
evaluate the integrands at the points at which data was taken, 


Weer) i1=1,2,...,N, with the subsequent application of 


7? 
Simpson's rule or the trapezoidal rule. Obviously, the inte- 
gration formulae used and the number of points at which the 
integrands are evaluated have an effect on the accuracies 


of the numerical integrals obtained. Nevertheless, it should 


be kept in mind that, no matter how accurate the integrals 
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are, the parameter estimates will be dependent on the co- 
efficients of the approximating polynomials used. Thus, too 
much accuracy in integration, which in turn means an increase 


in computing time, would not be justified. 


Ca Generation of Data 


1. Generation of Noise-Free Data 
For the generation of the numerical solution the infor- 
mation given in the previous pages was used. The remaining 


Operating variables are listed below: 


Zg = 0 

Ty = 1460°R 

D. = 4.026 in (inside diameter of tube) 

—d- = 5000 BTU/hr-sq-ft. (of the inside of tube) 
‘ 

P = 30 psia 

ines 60 Ibmoles/hr. 


With those values and initial conditions Equations 
(V-2) and (V-3) were integrated numerically, using step sizes 
of 10S 04 oiteand Of1 fo. hor each step 'saze neanly sidentical 
results were obtained, which were printed out at 25 ft. inter- 
vals. Smaller step sizes were not tried due to time limita- 
tions and possible round-off errors which might Sect 
ficant with the increasing number of calculations. The first 
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represent the noise-free measurements. These are presented 


in Table V-1. 


2. Generation of Noisy Data 


Twenty sets of noisy observations were obtained by the 
addition of random numbers with Gaussian distribution to the 


solution generated by numerical integration, using: 


~ * x 
Zi5 = 2; + ES Aas Tet Oh Oat ah ee 20 (V-22) 
= 7 1* Fes . ; 
Ti; api Tae di Ol ie ses. 90 4 EL) Cesc (V-23) 
* * 
where Zz. oT; = true observations of z and T at x = Xi 
aera? = noisy observations of z and T at x = x, 
J J fer tne-jth run: 


random numbers with Gaussian distribu- 
tion, zero mean, and standard deviations 
of O02 OF rand (5 20> vespectively:. 


As seen from Equation (V-22), the initial value of z 


was not subjected to noise. 


D.. Treatmenc oti vata 

Forsythe's method was used in fitting polynomials in 
x to all the observation sets. Comparison of variances 
corresponding to different polynomial degrees, led to the 
choice of fifth degree approximating polynomials for both z 


and T. In order to ensure that Zg = 0 -to an .accuyacy .of ,five 
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arg1tS; a weigntoon 10,000 was given ‘totz, ein: the fttting 


0 
procedure. Unit weights were assigned to the remaining 
points. In fitting polynomials to the temperature measure- 
ments, unit weights were assigned to all the 35 data points. 
The model errors and sensitivity coefficients were 
evaluated at the 35 original points at which data was taken, 


uSing the trapezoidal rule for integration. 


The initial parameter guesses were: 


The weighting matrix used was: 


Wo(NxN) 2(NxN) 


WONx2N) (V-24) 
Qonxn) “tr (Nxn) 
where W = ae (V-25) 
—z(NxN) Be —(NxN) 
Z 
W Bll (V-26) 
—T(NXN) we —(NXN) 
i 
LonxN) the identity matrix. (V-27) 


For the noisy observations oO, and o7 were set equal to 0.01 


and 5.0, respectively. For the noise-free observations, they 
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were set equal to 1.0 x 10 ° and 0.1, somewhat arbitrarily. 
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NOISE-FREE OBSERVATIONS 
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L. Results and Discussion 


1. Noise-Free Observations 
For convenience, the true values of the parameters 


are reproduced below: 


ay 397 .64¢ a, = 41.31 
Fifth, sixth, and seventh degree polynomials were tried 
as approximating functions. The results are presented in 
Table V-2. 
Studying Table V-2 shows that the estimates of a, are 
somewhat dependent on the degree of the approximating poly- 
nomial used, whereas the estimates of a, are relatively in- 


sensitive to the degree of the approximating polynomial. 


2. »Noisy Observations 

The results of the individual runs from one to twenty 
are shown in Table V-4. The summarized results are shown in 
Table V-3. 

Studying Tables V-3 and V-4, the following conclusions 
are drawn: 

(a) The estimates of a, are quite sensitive to observa- 
tion errors, whereas the estimates of a, are reasonably stable. 

(b) For each run, an increase in ay is accompanied by 


an increase in an. 


(c) The variance estimates of a, are quite sensitive 
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to the observation errors, and vary in proportion to ay (and 


vas 


a.) from run to run, whereas the variance estimates of as 


are quite stable. 


(d) The actual variances of the parameter estimates, 


Pe: ras 
S5 and ce » are larger than the averages of the individual 
1 Be 
parameter estimates, S¢ and $¢ 
a a. 


Conclusion (b) was directly supported by the high posi- 
tive values of the correlation coefficient p(a,,a,) observed 
during the individual least-squares runs. The correlation 


coefficient is defined by: 


p(a,.a,) 7 ooo Sane (V-28) 
a 
] 2 
where cov(a,,a,) is the estimate of the covariance of ay and 
an. 


It should be noted that the trends pointed out in 
Conclusions (a) to (c) emerge as a result of the characteris- 
tics of the system, i.e., the structure of the equations 
describing the system. Once the system is fixed by these 
equations, there is nothing that can be done to improve the 
results, no matter what kind of estimation technique is used, 
unless new experiments are run under more favorable condi- 
tions. Conclusion (d), on the other hand, can be attributed 
to the use of the approximating polynomials whose coefficients 


are treated as constants, rather than random variables. 


a 


bas) 8 ‘ohne 
<8 to 2ossmis2s sonst 


a 
" 1 
fF 


AAS 


ae ‘i ¢ 7 
ms Ad ¢? sited wre te hy oe te 


c2otemites retemsisq ant to. 
[subtvtbnt edd to 29pe78VvE A oe F . 
; _— 1 A - ea ny a 

' ; +8? ap 7 fo ¢ as : a 
-teoqg dptd sat yd bettoqgque vissevee 25W (4) an 
bsvysedo ( gb 84 gnstotttsoo ‘nots eT ST4es" $n ale 


norte feyi109 SAT .2nuy sSibbpesdbent” fsubtvibn 


xd bestteb ef nal 
i aS ry 7 ® Ads Pe 


bi ives Kf remepaaes >a then 
(88-¥) Pea = (ee ps) a 
ke crey q, ! « Van | 


5° 6 
ae | vl oe iz 


' brs 76 to SonsTYSBVOD sdd to stemrties sat af (gtr yadve9 aa a 
. say 

nt tuo bstntoq ebnsyt edt sent beton ad bfuode 31 | 

-2trsfo61edo oft to tfuesy 6 26 spteme (3) o+ (5) enotavfoned 

enofisups 9A to stusowrse sas ,.9.7 mat eye ond to°satek (nae : 

g2oit yd boxt? ef meteye oat 99nd mad eye edt ontdtyoeab | L 4 

esdy svo1qmt of snob 9d nea teat pntdson et stead | 2notseups itty 

.be2u ef suptaroad nottemites to bata teaw vottam on bg tUNane ry’, 
-fbnos sfdsvovst stom tabaw awn a8 iittala 

bajudtydts sd nebo ,baed isto eft ho .(b) ‘nofav fone? mest a. 

2instofttses szonw efstmonyloq pattemtxor1qgs pater rr “e. tT 

6 


2ofdstisy mobhen! asat 49Kde% ‘iene: ae 


ro. 


; Lat 
n't _ 
vm 


foes Ved Save : 
a : _ 


TABLE V-2 


PARAMETER ESTIMATES FROM NOISE-FREE OBSERVATIONS 


Degree of n x 
Approximating a, a. 
Polynomial 
5 48.18 41.12 
6 SO ae 483 41.36 
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TABLE V-3 


STATISTICAL RESULTS FROM NOISY OBSERVATIONS 
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TABLE V-4 


INDIVIDUAL ESTIMATES 


ay a> 4 
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F. Time Requirements 


Convergence was reached in 4 to 5 iterations per run 
on the average. The computing time required was 19 seconds 
per iteration. Polynomial fitting to the observations took 


about 13 seconds per data set. 


G. Comments on the Error Model 

With the error model discussed in previous pages on 
hand, one has various options in carrying out the integra- 
tions necessary to evaluate the model errors and sensitivity 
coefficients. As noted previously, the model errors and 
sensitivity coefficients were evaluated at 35 equally spaced 
points along the length of the reactor. Integration at a 
certain point was carried out making use of the integrand 
values at the previous points already calculated. By study- 
ing the profile of these integrands along the x axis (distance 
from the reactor inlet), one could decide to redistribute the 
points at which the integrands should be evaluated. Thus, 
more evaluations would be made in regions where the slopes 
of the integrands are relatively large. Still another 
possibility would be to evaluate the model errors and sensi- 
tivity coefficients at points other than those at which the 
observations were made, in addition to these. This would 
seem to be especially useful when the number of observations 
is small. However, one should not overlook the fact that, 
in this case, those values predicted by the approximating 


polynomial, rather than those actually observed, are used in 
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the evaluation of the model errors. Since the number of 
observations is small, one's confidence in the coefficients 
of the approximating polynomials would be relatively small, 
with the consequence that one's confidence in the values of 
the resulting integrals would be small, no matter how exact 
the integration scheme is. This limitation should always be 
kept in mind while estimating the reliability of the para- 


meter estimates. 


a 


| 


pe We i m 
| o* vodmun en 
anata Htta09 oda ae tno> 
iy we 
«! fame vovtiatey nie ane: | 


t96x9 nt veead on iene od fexypatnt pa 
9d 2eyeawls bluonde ROTSHS TET at amerae — 


ie a ey 


~s18q odd to yttitdstion ont iene att om 


pet reane we. : I ui i hy’! 
“a eh 


: sue e2 "ae 


P - y ae 
ise 7 e ‘ evaed? Te 
' ory 20% vara - 
| + Be 


: / sees Th gata 


‘ rT) ‘wd ae ie 


ady bi a 


i] 
‘i ia 
i 
‘a “>a a Rates: 
nieve 
IAB i 4 ae ae | 
oe | 
j ih vy | =) alo vege 
ed aS 


CHAPTER VI 


CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 


In this study, semi-empirical models have been deve- 
loped by introducing approximating functions into the integral 
transforms of the differential equations describing the 
system under study. In specific, ordinary power polynomials 
and exponential functions, expressed in terms of the independ- 
ent variable, have been considered. The former have been 
found to be more flexible in approximating a given set of 
observations. Moreover, they possess the desirable property 
that they can be used in approximating any kind of observa- 
tion set, regardless of the theoretical predictions about the 
shape of the profile of observations versus the independent 
variable. In some cases, however, it would be advantageous 
to use the knowledge about the theoretical predictions; if, 
say, theory predicts observations which will vary periodically 
with changes in the independent variable, then periodic func- 
tions of some kind should be considered as candidates for 
approximating functions. For cases where such trends are 
not expected, as was the case with the problems studied in 
this work, the use of ordinary power polynomials would prove 
of great advantage. 


Almost no numerical difficulties were encountered with 
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the initial parameter guesses picked prior to estimation. 
If the initial guesses do not lead to converging results, 
the data perturbation scheme proposed by Donnelly and Quon 
(9) can be incorporated into the estimation algorithm. 

The main goal of reducing the computation time due to 
integrations involved has been reached, not without a cost, 
however; the main problem now, is the difficulty encountered 
in establishing proper confidence intervals for the parameter 
estimates. Once the observations have been approximated by 
some functions, the parameters in these functions are con- 
Sidered to be constants for the rest of the computation. 

When the differential equations involved are linear, it is 
possible to express some parameters as functions of the para- 
meters which appear in the approximating functions. Subse- 
quently, confidence limits for the former can be established 
from a knowledge of the statistics of the latter. For non- 
linear differential equations, which represent more realistic 
situations, this kind of analysis is not possible. Therefore, 
the development of a general approach which would consider 
the parameters in the approximating functions as random 
variables with known statistics, and carry out the subsequent 
estimation accordingly, would be very desirable. 

In the approach developed in this work, systems consist- 
ing of only two simultaneous differential equations at the 
maximum have been considered. Another point of interest would 
be the investigation of the applicability of this approach to 


larger systems of differential equations. 
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NOMENCLATURE 


CHAPTER I 


p-vector of unknown parameters in the model 
function, defined by Equation (I-1) 

Nxp matrix of sensitivity coefficients 
function, defined by Equation (I-3) 

number of data points per run 

number of unknown parameters in the model 

pxp matrix, defined by Equation (1-4) 
objective function to be minimized, a function 


Of fel e 


p-vector, (2s Mails cise eee ss 
‘i 


independent variable and its initial value 
dependent variable and its initial value 
N-vector of observed values of y 


N-vector of predicted values of y by the 
algorithm 


p-vector of corrections to the current values 
of the parameters 
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N-dimensional vector of model errors (= y - y) 
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a scaler determining step length 


CHAPTER II 


p-vector of unknown parameters in the original 
model 


q-vectors of coefficients of the approximating 
function and their least-squares estimates 


Vector with dimension smaller than pt+q, ob- 
tained by combining the components of a and 

b in a manner determined by the specific pro- 
blem on hand 

Green's function for the initial value problem 
a function 

Green's function for the boundary value problem 
a function 

number of data points per run 


function, defined by Equation (II-5) 


number of unknown parameters in the original 
model 


number of coefficients of the AACE UT AS 
function 


constant appearing in the first boundary 
condition 


objective function to be minimized 


constant appearing in the second boundary 
condition 


a function 
a function 
independent variable 


initial and final values of x 
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dependent variable and its observed value at 
X=X. 
‘ 


function, defined by Equation (II-4) 


discrete and continuous model errors 


value of y'(x) at X=Xq 


constants appearing in the boundary conditions 


value of y(x) at X=Xo 
dummy variable, used for integration 


constants appearing in the boundary conditions 


function approximating the observations of y 
function obtained by substituting y(&) in the 


solution of the differential equation by the 
approximating function $(x:b) 


CHAPTER III 


mth coefficient of the polynomial approximat- 
ing c"(x) and its least-squares estimate 


Da-a,, and its least-squares estimate 
initialjguesssof De 


dimensionless measure of the composition of 
the reactant 


noise-free observation of c at X=X. 


noisy observation of c at X=X as for a certain 
runoand TOM shes athe run) in Specific 
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value of c at X=Xa for the jth run, obtained 


by evaluating the corresponding approximating 
polynomial of degree k 


Damkohler number and its least-squares estimate 


average of the individual least-squares esti- 
mates of Da 


number of runs with a certain data set 


cee degree of the polynomial approximating 
c( x 


number of data points per run 
reaction order and its least-squares estimate 


average of the individual least-squares esti- 
mates of n 


number of parameters estimated in Phase 2 


Peclet number and its least-squares estimate 


initial guess of Pe 


average of the individual least-squares esti- 
mates of Pe 


term accounting for the kinetics of decomposi- 
tion 


least-squares objective functions 


individual least-squares estimates of the 
variances of Da, Pe, and n 


variance estimates of Da, Pe, and n, obtained 
from their least-squares estimates Da, Pe, and n, 
USING. al lous 


averages of Save sae a using all runs 


measure of the performance of polynomial of 
degree k, in approximating the observations 
in nee en run 
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continuous and discrete weighting factors, 
reflecting one's confidence in each indivi- 
dual observation 


least-squares estimate of the mth coefficient 
of the Forsythe polynomial of degree m 


continuous and discrete model errors 


random number with Gaussian distribution and 
zero mean, used for generating noisy observa- 
tions 

dummy variable, used for integration 


standard deviation of the noise in the composi- 
tion measurements 


Forsythe polynomial of degree m 


CHAP TE Ray 


reactants 


functions of the unknown parameters, Ky» Kos 
and k, defined by Equation (IV-19) 


functions of the unknown parameters kK, >k,, and 
k defined by equation (IV-19) 


ig 

linear coefficients in Equation (IV-81) 
linear coefficients in Equation (IV-82) 
function, defined by Equation (IV-52) 


function, defined by Equation (IV-53) 


equilibrium constants defined by Equation 
(IV-11) 
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reaction rate constant for the decomposition 
ims pecies 1 to Species 4594 .J2lee oe Ito. 


eee 


reaction rate constants defined by Equation 
CIVeat2) 5 time”! 


k and k 


LoLtial guesses of Kk 9 3 


] + ] 
least-squares estimates of ky Kos and kK 3 


averages of the individual least-squares esti- 
mates of Ky > kK,» and k, 


optimum degree of polynomial P(t) 

optimum degree of polynomial Q(t) 

number of data points per run 

polynomial approximating the observations of xy 


mth coefficient of P(t) and its least-squares 
estimate 


polynomial approximating the observations of X 3 


mth coefficient of Q(t) and its least-squares 
estimate 


exponential function approximating the observa- 
tions of xy 


linear coefficients of R(t) 


exponential function approximating the observa- 
tions of X 3 


linear coefficients of S(t) 
least-squares objective function 


individual least-squares estimates of the 
variances of Ky» Ko > and kK 4 


variance estimates of Ky» Ko > and kK 35 obtained 


from their least squares estimates Ky> Kos and 
K35 Using ala uruns 


= 


not ssup3 vs bent Yeb sonssenog e3a% ROtIBAON 


ee 
go ne (agin aD ei er 
gi bas) XA Yo eodamives. 


-fi25 25%supe-tesal oe seaatom 


(+)9 fetmonyloq to sstpsb mumrsqo a 

(7)0 fetmonytoq hind i 

uy teq etntoa sisb Yo redmun 

rx to anottavisedo ott ontgsatiovens retmony tox 


2eo1supe-tesel 2tt bes (3)9 To snetotttsoo nim 
: atomises 


gx Yo enofssvisedo addi pntismixovggs Tetmonyloq 


2916upe-tesal att bas (2)0 to tastort¥soo Adam 
stemises 


“svigedo sds pritvemtxorggs nofrtoau fettnenogxs 
rx to enoht 


(3) to esnetaitteos bea dd 


-svisedo ond Abbie fut sottonut istduameuks 
oy de ghe TO) soge 

(3)2 ib 2tnstoatttsos ssantl 

erasiulYse—3 ane 


not donut avitostdo ze ssupe~-t 
at’. WHAT Tas jh 


ant to zotamttes 2evsupa-seeal Té 
ew brs al weed 
beantsido og bas « ep Yo. 264 5mt 


bNS «ot eA 2otemttes ae ae 


ee: a 
ime i 


% 


Me 


xy 9 My 3 


eo 220 o 20 


X 
Rees 


x 


<2 


Xx x 
dees 


sane 
less footed 


Bie oan 


124 


j os fied Sh 
2 3 


averages of § USminiGuea:): li KuNSs 


measure of time 


exponential coefficient of R(t) and its initial 
guess 


exponential coefficient of S(t) and its initial 
guess 


2N x 2N diagonal weighting matrix for the 
observations of xy and X. 


NXN diagonal weighting matrix for the observa- 
Clon OF Xy 


NxN diagonal weighting matrix for the observa- 
tion, OF X4 


normalized measures of the compositions of 
Remo and C 


values of Xa X and X9 at t=0 


o> 


equilibrium values of X, and X 4 


noise-free observations of xy and X 4 at t=t 


using observations of X and x 


at t=t., for a 
certain run 


é. 


using observations of X, and X 3 at t=t., for the 
Jen run. 7 spec. tae 


continuous model errors 
discrete model errors 


random numbers with Gaussian distribution and 
zero mean, used for generating noisy observa- 
tions 
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exponential coefficients in Equations (IV-81) 
and (IV-82) 


dummy variable, used for integration 


standard deviations of the noise in the obser- 
vations of X and X 4 


CHAPTER V 


reactant 
CrOss-seccional ared Of tube, Sd.it tc. 


unknown parameters 


initial guesses of a, and a, 


least-squares estimates of ay and a5 


averages of the individual least-squares esti- 
mates of ays and a. 


products 


heat capacities of the reactant and products, 
BTU/1lbmole-°R 


constant, defined by Equation (V-11) 


inside diameter of tube, ft. 
function, defined by Equation (V-17) 


' -] 
overall reaction rate term, sec 


linear and exponential coefficients in the 
overall reaction rate term 


optimum degree of polynomial P(x) 
optimum degree of polynomial R(x) 


number of data points per run 
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mean pressure inside the reactor, psia 

mth coefficient of P(x) 

heat input from furnace, BTU/hr-ft of tube 
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polynomial approximating the temperature 
observation 


mth coefficient of R(x) 
least-squares objective function 
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NxN diagonal weighting matrix for the tempera- 
ture observations 
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discrete model errors 
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APPENDIX 
A BRIEF OUTLINE OF THE NONLINEAR LEAST-SQUARES 


METHOD FOR PARAMETER ESTIMATION 


Given a model of the form: 
y = g(x:a) : CA-1) 


and a set of observations {Ys »X V=lo2c.8e,Nh, the estimation 


7? 
of the p-dimensional parameter vector a is desired. This is 


accomplished by minimizing: 
Laer 2 
oe! ) ee ) ly, ia g(x.:a)] (A-2) 


Starting with an initial set of parameter guesses, ag» 


the solution for the vector of corrections to the parameter 


guesses is given by: 


aes (A-3) 


where Go = Nxp matrix of sensitivity coefficients, 


= ¥ - Yo: 
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N-vector of predicted values of y by the algorithm 
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bo = p=vector of Corrections to the parameter guesses. 


In Equation (A-3), T stands for transpose. The improved set 


of parameters is obtained using: 
(A-4) 


The above procedure is repeated until every element of 6 is 
smaller than a predetermined error criterion. After this 
objective has been reached, the variance-covariance matrix 


Of parameters, Chase optained trom: 
C= of ee CAeay 


Bas ; f 
where q.. is the variance of the observation errors. In cases 
where an estimate of its value is not available, it can be 


approximated by: 


where a is the vector of the least-squares estimates of a, 
obtained with the algorithm outlined above. 

In cases where more than one dependent variable is 
observed, weights should be included in the analysis in order 
to express one's confidence in the observed values of the 


different variables. This is done in the following manner: 
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minere, this time G = Mxp matrix, of sensitivity coefficient, 
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M-vector residuals, 
M = NxK, 


K = number of dependent variables observed, 


W, 0 0 
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kK ae 49 929 eee 9 
K 
Lt = ONKN identity matrix, 
a = variance of the observation errors for the kth 


dependent variable. 


The diagonality of W follows from the assumption that observa- 


tions of different variables are uncorrelated. 
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